%% Portfolio Sorting: Quintiles (Table 4)
clear Q_ret_Spread_EW LS_ret_Spread_EW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INSERT SORTING VARIABLE AND RETURNS HERE
Sort_Var = (Stock_PRisk); 
returns = (call_deltahedge_returns_rb(1:end,:))*100;  
  % Note that data call and put option returns start from 02/2003, while 
  % PRisk start from 01/2003

% INSERT FACTORS FOR FACTOR MODELS
Factors1 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end)];
Factors2 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end) strdl_m(2:end) dvix_m(2:end)];
rf = rf_m(2:end)*100; % Risk-free rate
NW = 4; % Newey-West lag 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:size(Sort_Var,1)
    id1m = ~isnan(Sort_Var(i,:)+returns(i,:));  
    Returns = returns(i,:); Returns(:,find(id1m == 0)) = [];
    Betas_1M = Sort_Var(i,:); Betas_1M(:,find(id1m == 0)) = [];
      
    q_points_Spread_1M = prctile(Betas_1M, [0 20 40 60 80 100], 2);
    
    for j = 1: size(q_points_Spread_1M,2)-1
       Q_ret_Spread_EW(i,j) = nanmean(Returns(1,Betas_1M>q_points_Spread_1M(1,j)...
         & Betas_1M<=q_points_Spread_1M(1,j+1)));
       Q_Num_Obs(i,j) = sum(~isnan(Returns(1,Betas_1M>=q_points_Spread_1M(1,j)...
         & Betas_1M<=q_points_Spread_1M(1,j+1))));
    end   
end

Obs = nanmean(Q_Num_Obs);

Result_SingleSort = nan(1,16);
% Equal-weighted returns
Result_SingleSort(1,1:2:9) = nanmean(Q_ret_Spread_EW);
LS_ret_Spread_EW = Q_ret_Spread_EW(:,5)-Q_ret_Spread_EW(:,1); 
Result_SingleSort(1,11) = nanmean(LS_ret_Spread_EW);
for j = 1:5
    stats_Spread_EW = nwest(Q_ret_Spread_EW(:,j)-rf,ones(size(Q_ret_Spread_EW(:,j))),NW);
    Result_SingleSort(1,j*2) = stats_Spread_EW.tstat;
end
stats_Spread_EW = nwest(LS_ret_Spread_EW,ones(size(LS_ret_Spread_EW)),NW); 
Result_SingleSort(1,12) = stats_Spread_EW.tstat;

stats_Spread_EW = nwest(LS_ret_Spread_EW,[ones(size(LS_ret_Spread_EW)) Factors1],NW); 
Result_SingleSort(1,13) = stats_Spread_EW.beta(1);
Result_SingleSort(1,14) = stats_Spread_EW.tstat(1);

stats_Spread_EW = nwest(LS_ret_Spread_EW,[ones(size(LS_ret_Spread_EW)) Factors2],NW); 
Result_SingleSort(1,15) = stats_Spread_EW.beta(1);
Result_SingleSort(1,16) = stats_Spread_EW.tstat(1);

Result_SingleSort = Result_SingleSort';

clear Sort_Var returns Factors rf NW i j k stats_Spread_EW Q_ret_Spread_EW returns ...
      stats_Spread_FFC_EW stats_Spread_VW Factors Betas_1M q_points_Spread_1M ...
      LS_ret_Spread_EW Factors1 Factors2 Q_Num_Obs Returns

%% Quintile Analysis - OptionValue Weighted
clear Q_ret_Spread_VW LS_ret_Spread_VW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INSERT SORTING VARIABLE AND RETURNS HERE
Sort_Var = (Stock_PRisk); 
returns = (call_deltahedge_returns_rb(1:end,:))*100; 
Weight = MVO; % Option market value: Open interest * price
  % Note that data call and put option returns start from 02/2003, while 
  % PRisk start from 01/2003

% INSERT FACTORS FOR FACTOR MODELS
Factors1 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end)];
Factors2 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end) strdl_m(2:end) dvix_m(2:end)];
rf = rf_m(2:end)*100; % Risk-free rate
NW = 4; % Newey-West lag 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for i = 1:size(Sort_Var,1)
    id1m = ~isnan(Sort_Var(i,:)+returns(i,:)+Weight(i,:));  
    Returns = returns(i,:); Returns(:,find(id1m == 0)) = [];
    W = Weight(i,:); W(:, find(id1m == 0)) = [];
    Betas_1M = Sort_Var(i,:); Betas_1M(:,find(id1m == 0)) = [];
      
    q_points_Spread_1M = prctile(Betas_1M, [0 20 40 60 80 100], 2);
    
    for j = 1: size(q_points_Spread_1M,2)-1 
       Q_ret_Spread_VW(i,j) = Returns (1,Betas_1M > q_points_Spread_1M(1,j) & Betas_1M...
         <= q_points_Spread_1M(1,j+1))*(W(1,Betas_1M > q_points_Spread_1M(1,j)...
         & Betas_1M <= q_points_Spread_1M(1,j+1))./nansum(W(1,Betas_1M >...
         q_points_Spread_1M(1,j) & Betas_1M <= q_points_Spread_1M(1,j+1))))';
    end   
end

Result_SingleSort = nan(1,16);
% For value-weighted
Result_SingleSort(1,1:2:9) = nanmean(Q_ret_Spread_VW);
LS_ret_Spread_VW = Q_ret_Spread_VW(:,5)-Q_ret_Spread_VW(:,1); 
Result_SingleSort(1,11) = nanmean(LS_ret_Spread_VW);
for j = 1:5
    stats_Spread_VW = nwest(Q_ret_Spread_VW(:,j),ones(size(Q_ret_Spread_VW(:,j))),NW);
    Result_SingleSort(1,2*j) = stats_Spread_VW.tstat;
end
stats_Spread_VW = nwest(LS_ret_Spread_VW,ones(size(LS_ret_Spread_VW)),NW); 
Result_SingleSort(1,12) = stats_Spread_VW.tstat;

stats_Spread_VW = nwest(LS_ret_Spread_VW,[ones(size(LS_ret_Spread_VW)) Factors1],NW);
Result_SingleSort(1,13) = stats_Spread_VW.beta(1);
Result_SingleSort(1,14) = stats_Spread_VW.tstat(1);

stats_Spread_VW = nwest(LS_ret_Spread_VW,[ones(size(LS_ret_Spread_VW)) Factors2],NW);
Result_SingleSort(1,15) = stats_Spread_VW.beta(1);
Result_SingleSort(1,16) = stats_Spread_VW.tstat(1);

Result_SingleSort = Result_SingleSort';

clear id1m id1M id1MM Returns W Betas_1M q_points_Spread_1M Sort_Var ...
      i j k stats_Spread_EW Q_ret_Spread_W Q_ret_Spread_W returns Weight...
      stats_Spread_FFC_VW stats_Spread_FFC_EW stats_Spread_VW Factors rf NW ...
      LS_ret_Spread_W LS_ret_Spread_W Factors1 Factors2 rf

%% Quintile Analysis - capped weighted
clear Q_ret_Spread_VW LS_ret_Spread_VW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INSERT SORTING VARIABLE AND RETURNS HERE
Sort_Var = (Stock_PRisk); 
returns = (call_deltahedge_returns_rb(1:end,:))*100; 
Weight = Stock_CAP_capped; 
  % Note that data call and put option returns start from 02/2003, while 
  % PRisk start from 01/2003

% INSERT FACTORS FOR FACTOR MODELS
Factors1 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end)];
Factors2 = [mktrf_m(2:end) smb_m(2:end) hml_m(2:end) umd_m(2:end) strdl_m(2:end) dvix_m(2:end)];
rf = rf_m(2:end)*100; % Risk-free rate
NW = 4; % Newey-West lag 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% Long-short strategy: Value weighted and Equally weighted (sorted by Spread)
for i = 1:size(Sort_Var,1)
    id1m = ~isnan(Sort_Var(i,:)+returns(i,:)+Weight(i,:));  
    Returns = returns(i,:); Returns(:,find(id1m == 0)) = [];
    W = Weight(i,:);W(:, find(id1m == 0)) = [];
    Betas_1M = Sort_Var(i,:); Betas_1M(:,find(id1m == 0)) = [];
      
    q_points_Spread_1M = prctile(Betas_1M, [0 20 40 60 80 100], 2);
    
    for j = 1: size(q_points_Spread_1M,2)-1 
       Q_ret_Spread_VW(i,j) = Returns (1,Betas_1M > q_points_Spread_1M(1,j) & Betas_1M...
         <= q_points_Spread_1M(1,j+1))*(W(1,Betas_1M > q_points_Spread_1M(1,j)...
         & Betas_1M <= q_points_Spread_1M(1,j+1))./nansum(W(1,Betas_1M >...
         q_points_Spread_1M(1,j) & Betas_1M <= q_points_Spread_1M(1,j+1))))';
    end   
end

Result_SingleSort = nan(1,16);
% For value-weighted
Result_SingleSort(1,1:2:9) = nanmean(Q_ret_Spread_VW);
LS_ret_Spread_VW = Q_ret_Spread_VW(:,5)-Q_ret_Spread_VW(:,1); 
Result_SingleSort(1,11) = nanmean(LS_ret_Spread_VW);
for j = 1:5
    stats_Spread_VW = nwest(Q_ret_Spread_VW(:,j),ones(size(Q_ret_Spread_VW(:,j))),NW);
    Result_SingleSort(1,2*j) = stats_Spread_VW.tstat;
end
stats_Spread_VW = nwest(LS_ret_Spread_VW,ones(size(LS_ret_Spread_VW)),NW); 
Result_SingleSort(1,12) = stats_Spread_VW.tstat;

stats_Spread_VW = nwest(LS_ret_Spread_VW,[ones(size(LS_ret_Spread_VW)) Factors1],NW);
Result_SingleSort(1,13) = stats_Spread_VW.beta(1);
Result_SingleSort(1,14) = stats_Spread_VW.tstat(1);

stats_Spread_VW = nwest(LS_ret_Spread_VW,[ones(size(LS_ret_Spread_VW)) Factors2],NW);
Result_SingleSort(1,15) = stats_Spread_VW.beta(1);
Result_SingleSort(1,16) = stats_Spread_VW.tstat(1);

Result_SingleSort = Result_SingleSort';

clear id1m id1M id1MM Returns W Betas_1M q_points_Spread_1M Sort_Var ...
      i j k stats_Spread_EW Q_ret_Spread_W Q_ret_Spread_W returns Weight ...
      stats_Spread_FFC_VW stats_Spread_FFC_EW stats_Spread_VW Factors rf NW ...
      LS_ret_Spread_W LS_ret_Spread_W Factors1 Factors2 rf
